function [rdat,xcoord,ycoord] = rscan(M0,varargin)

%%% Computes a radial scan for use in automated computation of macular pigment maps.
%%% Modified by Emily King, based on the function found at
%%% http://www.mathworks.com/matlabcentral/fileexchange/18102, created by
%%% Narupon Chattrapiban
%%% Last modified: 27 July 2010

%%% Some of his documentation:
% RDAT = RSCAN(M0,VARARGIN)
% Get radial scan of a matrix using the following procedure:
% [1] Get coordinates of a circle around an origin.
% [2] Average values of points where the circle passes through.
% [3] Change radius of the circle and repeat [1] until rprofile is obtained.
%
% Draw Circle: 
% [ref] http://www.mathworks.com/matlabcentral/fileexchange/
%       loadFile.do?objectId=2876&objectType=file

%%% My documentation:
% rlim = max radius, default maximum possible
% xavg = center column, default center of image
% yavg = center row, default center of image
% dispflag = 0/1 determines whether original image displayed (This is
% slightly different than the original function, where this controled both
% displays.), default on
% [Removed from original function] dispflagc = 0/1 determined whether this
% stupid red ball would trace all of the circles in the radial profile.
% Totally pointless and slow.
% [Not in original function] dispflag2 = 0/1 determines whether the plot of
% the radial profile is displayed, default on
% squeezx = a (multiplicative) deformation of x to obtain an elliptical scan, default 1
% squeezy = a (multiplicative) deformation of y to obtain an elliptical
% scan, default 1
% rot = a rotation amount in radians (useful for elliptical scans), default
% 0
% rstep = the steps taken in radius size, default 1

% rdat = sequence of the average values over each circle of growing radius
% [xcoord ycoord] = actual pixels of the final circle

% Example usage:
% rscan(Image);
% [rdata, col, row] = rscan(Image,'xavg',100,'dispflag',0,'squeezy',0.8);




if nargin < 1, 
    disp('This is a Demo');
    M0 = 'demo';
end

if strmatch('demo',lower(M0),'exact')
    %M0 = peaks(200);
    [xx,yy] = meshgrid(linspace(-3,3,201));
    [phi,rho] = cart2pol(xx,yy);
    M0 = besselj(1,rho);
    clear xx yy phi rho;
end

xavg = size(M0,2)/2;
yavg = size(M0,1)/2;
Rlim = floor(min(size(M0))/2)-1;
dispFlag = 1;
dispFlag2 = 1;
%dispFlagC = 0;
rot = 0; %% radian
squeezx = 1; %% 0.80;
squeezy = 1;
rstep = 1;

if exist('varargin','var')
    L = length(varargin);
    if rem(L,2) ~= 0, error('Parameters/Values must come in pairs.'); end
    for ni = 1:2:L
        switch lower(varargin{ni})
            case 'xavg', xavg = varargin{ni+1};
            case 'yavg', yavg = varargin{ni+1};
            case 'squeezx', squeezx = varargin{ni+1};
            case 'squeezy', squeezy = varargin{ni+1};
            case 'rot', rot = varargin{ni+1};
            case 'rlim', Rlim = varargin{ni+1};
            case 'dispflag', dispFlag = varargin{ni+1};
            %case 'dispflagc', dispFlagC = varargin{ni+1};
            case 'dispflag2', dispFlag2 = varargin{ni+1};
            case 'rstep', rstep = varargin{ni+1};
            otherwise error(strcat('''',varargin{ni},'''',' is an unknown argument type.'));
        end
    end
end


yxz = size(M0);
Rbnd = floor(min(yxz)/2)-1;
if Rlim > Rbnd, Rlim = Rbnd; end
    
for nRho = 1:rstep:floor(Rlim),
NOP = round(2*pi*nRho);
THETA=linspace(0,2*pi,NOP);
RHO=ones(1,NOP)*round(nRho);
[X,Y] = pol2cart(THETA,RHO);
X = squeezx*X;
Y = squeezy*Y;
[THETA,RHO] = cart2pol(X,Y);
[X,Y] = pol2cart(THETA+rot,RHO);
X = X + xavg;
Y = Y + yavg;

if dispFlag,
    h1 = figure(100);clf;box on;
    set(h1,'position',[10 500 400 300]);
    set(h1,'units','degrees'); 
    set(gca,'units','degress'); 
    imagesc(M0);axis image;hold on;
    % H = plot(X,Y,'c-');
    line([xavg xavg],[1 yxz(1)],'color',[1 1 0],'linewidth',1); 
    line([1 yxz(2)],[yavg yavg],'color',[1 1 0],'linewidth',1);
end

%%%
X = round(X);
Y = round(Y);
%%%

clear dat uxy pxy mxy mx nx my ny;
dat = [X;Y];
uxy = diff(dat,1,2);
uxy = [[1;1],uxy];
pxy = union(find(uxy(1,:)~=0),find(uxy(2,:)~=0));
dat = dat(:,pxy);
rdat(nRho) = sum(M0(yxz(1)*(dat(1,:)-1) + dat(2,:)))/length(dat);

if dispFlag,
    H = plot(dat(1,:),dat(2,:),'y-');
%    if dispFlagC,
%        for nrn = 1:length(dat)
%            H = plot(dat(1,nrn),dat(2,nrn),'m.','MarkerSize',12);
%            drawnow;
%            delete(H);
%        end
%    end
    drawnow;
end
end

% rdat = rdat/max(rdat);  Removed from original function because we want
% quantitative values.
xcoord = dat(1,:);
ycoord = dat(2,:);

if dispFlag2,
% VNR(undo) h2 = figure(102); clf; hold on;
% VNR(undo) set(h2,'position',[10 100 400 300]);
% VNR(undo) plot(rdat,'rx','linewidth',2);axis tight;
% plot(M0(round(xavg)+1,round(yavg)+1:end),'b','linewidth',2);
end



